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SscJon  I.  INTRODUCTION 


This  report  documents  the  results  of  a study  on  inertial  measurement 
unit  (DflJ)  self-alignment  techniques  for  PERSHING  II.  The  primary  concern 
is  the  fine  alignment  of  a fixed  base  inertial  platform  where  the  base  is 
subjected  to  ground  vibration  and  wind  buffeting. 

A thorough  discussion  will  be  made  on  the  relationships  among  drifts 
and  misalignments  of  a platform.  A gyrocompassing  equation  will  be 
derived.  S?.  serai  concepts  useful  for  forming  self-alignment  procedures 
will  be  isc*vised  vlth  the  help  of  developed  analytics. 

A new  lea***  square  regression  algorithm,  specially  for  IMU  alignment, 
"'ll  be  develo  1.  The  superiority  of  this  algorithm  will  be  demonstrated 
.rough  theoretical  analysis,  experimental  results,  and  hypothetical 
examples. 

The  scope  of  this  study  can  best  be  seen  from  the.  Table  of  Contents. 
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Section  II.  IMU  SELF-ALIGNMENT 


An  IMU  can  be  aligned  Co  an  earth  fixed  coordinate  system  at  most 
latitudes  on  the  earth's  surface  by  using  the  information  derived  from 
the  output  of  the  unit's  sensors.  This  type  alignment  of  an  IMU  is  cal- 
led leveling  and  gyrocompassing.  Two  fundamentally  different  approaches 
are  used  to  accomplish  this:  (1)  the  gimbaled  platform  of  the  IMU  is 

physically  driven  to  align  with  the  earth  coordinates , and  (2)  the  align- 
ment is  achieved  analytically  by  determining  the  misalignments  of  platform 
axes  with  respect  to  the  earth  coordinates.  The  second  approach  has  the 
advantage  of  faster  gyrocompassing,  but  at  the  expense  of  a high  speed 
digital  computer.  Our  present  study  is  centered  on  the  second  approach, 
namely,  the  "analytic  gyrocompassing". 

The  earth  fixed  coordinate  system  adopted  in  this  study  is  shown  in 
Figure  1,  There  the  three  orthogonal  axes  are  N,  E,  and  A representing 
north,  east,  and  azimuth,  respectively.  As  a result  of  this  choice  of 
coordinates,  the  azimuth  component  of  the  earth  rate  is  a negative  quan- 
tity as  shown  in  Figure  2. 


Figure  1.  Earth  fixed  coordinates. 


Figure  2.  Earth  rate  components. 


Figure  3 depicts  the  involvement  of  the  analytic  gyrocompassing.  A 
large  amount  of  data,  obtained  from  the  outputs  of  platform  sen* res,  is 
reduced  to  a set  of  parameters  through  a chosen  data  reduction  technique 
Then,  the  platform’s  azimuth  is  determined  from  these  parameters  using  a 
gyro compassing  equation. 


It  is  obvious  that  better  hardware  allows  more  accurate  gyrocompas- 
sing. For  a given  hardware  system,  different  gyrocompassing  procedures 
can  be  formed  by  different  combinations  of  basic  concepts.  Thus  the 
accuracy  of  gyrocompassing  depends  on  several  factors: 


1)  Platform  hardware 

2)  Gyrocompassing  procedure 


3)  Gyrocompassing  equation 

4)  Data  reduction  algorithm 

5)  Computer  dependent  errors. 
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Figure  3.  Analytic  gyrocompassing. 


Figure  4 shows  a general  model  of  a platform  with  error  sources 
indicated.  Each  solid  line  represents  a hardware  connection  while  each 
dashed  line  represents  the  path  of  a sensed  signal. 


The  •'electronic  and  network"  block  can  be  implemented  for  various 
purposes  such  as  maintaining  the  platform  at  a specific  orientation, 
compensating  for  errors,  and  improving  the  platform  dynamics.  Our  pre- 
sent purpose  is  to  align  the  platform  coordinates  to  the  earth  fixed 
coordinates. 
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Section  Ilf.  PLAT  FOR  W DRIFTS  [1] 


The  analytic  model  needed  for  data  reduction  and  the  gyrocompassing 
equation  needed  for  azimuth  determination  are  derived  from  the  drift 
characteristics  of  the  platform.  Therefore  understanding  the  drift  char- 
acteristics is  a prerequisite  to  the  development  of  gyrocompassing 
techniques . 


For  a gimbaled  platform,  the  platform  axes  are  slaved  to  the  gyros. 
The  time  constants  of  platform  servos  are  usually  much  smaller  than  the 
gyrocompassing  time  (on  the  order  of  milliseconds  versus  minutes).  Hence 
the  gyro  drift  contributes  instantaneously  to  the  platform  drift  of  the 
same  amount.  Thus  the  terms  "platform  drift"  and  "gyro  torquing  rate" 
become  synonymous. 


1 . Kinematics  of  Platform  Drift 

Consider  a platform  which  has  been  coarsely  aligned  to  the 
earth  fixed  coordinates.  The  deterministic  torquing  rate  for  each  gyro 
consists  of  the  self-axis  earth  rate  component,  Che  cross-axis  earth  rate 
component,  and  the  gyro  drift. 

Let  6„  - misalignment  about  north  axis 

S* 

0 = misalignment  about  east  axis 

L 

6 ^ - misalignment  about  azimuth  axis 
= north  gyro  drift 
D_  = east  gyro  drift 

£i 

«*  azimuth  gyro  drift  race 

ft  * earth  rate 

L * latitude  of  launchsite 

ft^  = ft  cos  L = north  component  of  earth  rate 

ft^  = -ft  sin  L = azimuth  component  of  earth  rate 


Kjj  = north  gyro  torquer  scale  factor  error 
” azimuth  gyro  torquer  scale  factor  error. 


-~3i 
. -A 


8 


-^fe-Vaq&A^* 


feSST- 


The  torquing  rate  for  each  gyro  is  obtained  as  follows: 
for  the  north  gyro, 

§N  " % + °N  " flA  8in  eE  " V1  " C°S  6 A COS  V + h % (1> 

where  0.  sin  0 = cross-axis  earth  rate  due  to  misalignment 

A Ci 

f>  (1  - cos  0.  cos  0 ) = change  of  self-axis  earth  rate  due  to 

misalignment 

^ 53  rate  due  to  torquer  scale  factor  error; 
for  the  east  gyro, 

«E  ' °E  + "A  8ln  SN  - °»  Sin  «A  (2) 

where  sin  0^  - sin  0^  = cross-axis  earth  rates  due  to  misalignment; 

for  the  azimuth  gyro, 

" nA  + D*  + “»  Sl"  SE  - V1  - COS  SH  COS  V + KA  SiA  <3) 

where  sin  0£  * cross-axis  earth  rate  due  to  misalignment 

n.(l  - cos  0 cos  0 ) * change  of  self-axis  earth  rate  due  to 

misalignment 

= rate  due  to  torquer  scale  factor  error. 

From  Equations  (1),  (2),  and  (3)  the  non-nominal  parts  of  the 
torquing  rates  for  north,  east,  and  azimuth  gyros  (or,  equivalently,  the 
drifts  for  north,  east,  and  azimuth  axes  of  the  platform)  are, 
respectively, 

- °N  - “a  8in  eE  ’ V1  - CCS  COS  6E>  + h%  (4) 

°E  = °E  + QA  8in  0N  ' ‘7N  8in  6 A (5) 

°A  = °A  + aN  Sin  eE  ’ aA(1  ’ C0S  eK  C°8  ^E^  + KA  aA  * (6) 

If  the  misalignment  0 , 0 , and  0.  are  sufficiently  small,  small 
angle  approximations  for  sine  and  cosine  functions  can  be  used.  Under 


(5) , and  (6)  reduce  to 

-°AeE+*MSi 

(7) 

-°N0A+ftA0N 

(8) 

+ 0E  + **  ^ * 

(9) 

E E 


Monitoring  values  of  calibrated  gyro  torquing  currents  provide  a way  of 
determining  the  platform  drifts. 


2.  Time  Functions  of  Drifts 

Equations  (7) , (8) , and  (9)  show  that  the  drifts  along  three 
platform  axes  are  coupled  together  by  the  misalignments  0^,  0g,  and  0 

Understanding  this  coupling  effect  is  important  to  accurate  gyrocompassing. 

It  is  reasonable  to  assume  that  during  the  period  of  gyrocompassing, 

0.,  the  azimuth  misalignment  is  constant.  With  this  in  mind.  Equations 
A 

(7)  and  (8)  can  be  further  developed  into 


°N  " ^A^EO  + / DECT)dT^  + ^ ^ 

t 

DN-nA0EO+KNS  *nA  / DE<T)dT 


DE  “ °E  ’ 0A  + fiA  0N 


°E  ’ °N  0A  + nA<0NO  + / DN(T)dT 

0 

t 

°E  “ S 0A  + nA  eN0  + aA  / DN(T)dT 


By  defining 


deo  - de 


Si  0A  + nA  0NO  “ DE  " Si  eA 


DN0  “ DN  * nA  0EO  + *11  Si 


(12) 


which  represent  the  initial  values  of  the  drift  for  east  and  north  axes, 
the  two  drift  equations  become 


Di  + ftA  J DE<T>dr  * °N0 


DE-«A  / DN*T*dT  “ °E0  * 


(13) 


(14) 


The  time  functions  D*  (t)  and  D* (t)  can  be  solved  from  Equations  (13) 

fi  £ 

and  (14)  using  transform  method.  Taking  the  Laplace  transform  of  both 
equations. 


>;<»>  + t de<» 


HO 


de<8>  - T Di<8>  - ¥ 


(15) 

(16) 


Solving  fqj:  D^(s)  and  D^(s)  yields 


0^(a)  - D, 


8 « 

NO  _2  , J " °E0  2 


s + ft. 


D£(s)  = D£0  2 


-£~^+  D,’ 


s + ft. 


ft. 


S + ft. 


N0  *2  + »l 


(17) 


(18) 


By  taking  the  inverse  Laplace  transform  of  Equations  (17)  and  (18),  the 
corresponding  time  functions  are,  respectively. 


DN(t)  “ °N0  008  V " DE0  8in  V 

DE(t)  * DE0  008  nAfc  + °N0  8ln  flAfc 


(19) 

(20) 


I 
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In  matrix  form, 


d;<« 


cos  fi^t  - sin 


sin  nAt  cos 


which  shows  that  the  vector  drift  at  any  time  is  a rotation  of  angle  fit 
from  its  initial  vector  drift. 


For  small  value  of  ft. t, 

A 


„ 2 2 

cos  Sl^t  = 1 - — — 


sin  nAt  =s  nAt 


Then  Equations  (17)  and  (18)  can  be  approximated  by 

2 2 

"i<*>  ■ Dso  - deo  V - dmo  V-  <2 

n2  t2 

DE(t>  - DE0  + “ho  V - DE0  -V  • 

If  the  second  order  terms  are  negl igible,  Equations  (23)  and  (24)  can 
further  be  approximated  by 

DN(t)  = dn0  - d£0  nAt  (2 

DE(t>  " °E0  + DN0  V * <2 


3.  Time  Functions  of  Misalignments 

With  the  help  of  Equations  (19)  and  (20),  the  misalignments 
0„(t)  and  0_(t)  can  be  determined  as  follows: 


Section  IV.  GYROCOMPASSING  EQUATION 


Gyro compa s s ing , the  determination  of  azimuth  misalignment*  requires 
the  knowledge  of  drifts  and  misalignments  along  the  north  and  east  plat- 
form axes.  A gyrocompassing  equation  and  its  accuracy  are  discussed  here, 


1.  The  Equation 

By  solving  Equation  (8)  for  0 . (t) , the  azimuth  misalignment  at 
any  time  t is  obtained  as 


eA<t> 


PE  - D'(t)  + aA  yt) 

“’ll 


(29) 


Substituting  the  details  of  D^(t)  and  from  Equations  (18)  and  (25) 

into  Equation  (29),  all  sine  and  cosine  terms  cancel.  The  result  is  the 
"gyrocompassing  equation"  sought. 


de  • 

v-«  - — 


EO 


flA  °N0 


(30) 


Intuitively,  it  can  also  be  said  that  Equation  (30)  comes  directly 
from  Equation  (29)  since 


SA<t>  - 6A0  = 


°E  " °E0  + nA  °N0 


“n 


2.  Gyrocompassing  Accuracy 

The  ultimate  gyrocompassing  accuracy  is  limited  by  the  follow- 
ing uncertainties: 

ADg  = East  gyro  drift  uncertainty 

ACg  * Uncertainty  in  platform  drift  about  its  east  axis 

* Uncertainty  in  the  initial  platform  misalignment  about 
its  north  axis. 

From  Equation  (30)  the  uncertainty  in  gyrocompassing  is  obtained  as 

adp  - adJ  + n. 


iipp 


Since  the  sign  of  ea'ti  individual  uncertainty  is  not  known,  an  upper 
bound  of  the  gyrocoinpassing  error  is  given  by 


Itf.l  s 


iADg  - AD.  I +'|a»AeM„| 


E 


A NO1 


°N 


(32) 


In  Equation  (32) , A and  A^q  are  in  radians  while  ADg,  ADg,  0^,  and  5^ 
have  the  same  unit.  I 
should  be  replaced  by 


have  the  same  unit.  If  A 0A  and  a©n0  are  in  arcseconds.  Equation  (32) 


| AD'  - ADp| 

|a0a|  5 206,280  ^ S-+  lASyo  tan  L| 


(33) 


“A Si 

where  L is  the  launchsite  latitude. 

Consider  an  example  where 

L » 45  degrees 

AD£  - ADg  * 0.003  deg/hr 

A0„„  * 2 arcsec 

NU 

Since  0 » 15  degrees/hour, 

Ojj  ■ 0 cos  L ■ 10.61 

tan  L = 1 . 

Equation  (33)  gives 

|a©aI  S 1 206,280  X | + 2 = 58.3  + 2 = 60.3  arcsec 

Notice  that,  in  this  example,  the  platform  east  axis  drift  uncertainty 
contributes  most  of  the  error  in  gyrocompassing. 
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Section  V.  CONCEPTS  FOR  SELF  ALIGNMENT 


This  section  presents  a discussion  or  several  useful  concepts 
which  can  be  chosen  to  form  different  IMU  self-alignment  procedures. 


1 . Gyro  Drift  Determination 

Determination  of  gyro  drifts  using  the  information  within  the 
atfonu  system  is  also  called  "autobiasing. " Here  we  shall  be  con- 
.ned  with  drifts  about  north  and  east  axes.  There  are  two  different 
methods  of  autobiasing  gyros,  namely,  the  "closed-loop  method"  and  the 
"open-loop  method."  The  choice  between  the  two  depends  on  the  relative 
uncertainty  between  the  gyro  torquer  scale  factor  error  and  the  acceler- 
ometer scale  factor  error. 

Referring  to  Figure  4,  the  platform  alignment  system  consists  of 
two  Schuler  loops.  For  the  closed-loop  method,  both  Schuler  loops  are 
closed  and  gyros  are  torqued  at  the  rates  given  by  Equations  (1),  (2), 
and  (3).  Under  the  fine  alignment  condition,  platform  is  sufficiently 
level  such  that  small  angle  approximations  for  trigonometric  functions 
are  satisfactory.  Therefore  Equations  (7)  and  (8)  give  the  non-nominal 
torquing  rate  for  the  north  and  east  gyros.  In  general,  there  are 
biases  in  accelerometers,  so  the  terms  and  0^0^  may  not  be  small. 

However,  In  all  practical  cases,  the  biases  are  known.  Therefore 
their  effect  on  platform  drift  is  known.  Thus  it  can  be  said  that  the 
difference  between  and  the  corresponding  accelerometer  bias  effect 

is  small,  and  between  and  its  corresponding  accelerometer  bias 

effect  is  also  small.  Under  this  condition.  Equations  (7)  and  (8)  are 
further  reduced  to 


°N  ~ °N  + KNnN 


°E  “ °E  * Va 


(34) 

(35) 


The  quantities  and  are  obtained  by  measuring  the  torquing 

currents  of  north  and  east  gyros.  Assuming  that  is  known,  the  north 

gyro  drift  can  be  accurately  determined  if  K^,  the  torquer  scale 

factor  error,  is  known.  However,  the  uncertainty  in  (1^6 A is,  in  general, 

so  large  that  there  is  no  way  to  accurately  determine  the  east  gyro 
drift  D£  from  D^. 

Often  the  knowledge  of  Kjj  is  not  available  to  the  degree  of  pre- 
cision desired.  Under  this  condition,  accurate  and  rapid  determination 
of  D^  from  the  closed-loop  information  is  difficult. 
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The  effect  of  uncertainty  in  north  gyro  torquer  scale  factor  error 
can  be  eliminated  entirely  by  not  torquing  the  north  axis  of  the  plat- 
form physically.  Instead,  an  analytically  torqued  north  axis  is  main- 
tained in  the  computer  by  on-line  computation.  This  method  is  called 
"zero- torquing  measurement."  Zero-torquing  is  accomplished  by  opening 
Schuler  loops  at  places  Indicated  in  Figure  4.  Therefore  the  method 
is  an  "open-loop  method."  Under  this  condition,  platform  level  is  not 
maintained,  so  accelerometers  receive  larger  inputs.  Thus  the  uncer- 
tainty in  the  accelerometer's  scale  factor  error  becomes  more  important. 


In  the  open-loop  method,  measurements  are  taken  at  outputs  of  both 


accelerometers.  From  the  measurements,  and  are  determined. 


Because  of  zero-torquing,  plays  no  part  in  drift  determination,  so 
Equation  (34)  becomes 


dk  “ dn  • 


(36> 


which  is  an  attractive  way  to  determine  D^.  However,  is  still  given 
by  Equation  (35)  where  separation  of  Dg  from  -nN©A  is  difficult. 


To  conclude,  it  is  seen  that  whether  the  closed-loop  method  or 
open-loop  method  is  used  to  determine  gyro  drifts,  only  the  north  gyro 
drift  can  be  accurately  determined.  Reference  1 contains  several 
numerical  examples  to  illustrate  this  phenomenon.  The  technique  of 


determining  D'  and  D'  from  the  measurements  for  the  open-loop  method 
N E 


will  be  discussed  in  detail  in  Sections  VI  and  VII. 


Two-Position  Gyrocompassing 


Recall  from  Equation  (30)  that  an  accurate  gyrocompassing  can 
be  achieved  only  if  an  accurate  determination  of  the  east  gyro  drift  D£ 

can  be  obtained.  Since  accurate  drift  determination  can  be  made  only 
for  north  gyro,  a two-position  scheme  can  be  devised  to  take  this  advan- 
tage. The  scheme  may  consist  of  the  following  steps: 


a)  Autobiasing  the  north  gyro  to  determine  D 


N 


b)  Slewing  the  platform  approximately  90  degrees  so  that 
the  north  gyro  becomes  an  east  gyro,  and  the  original  east  gyro  becomes 
a south  gyro 

c)  Gyrocompassing  is  performed  with  the  present  east  gyro 
whose  drift  has  been  accurately  determined,  enabling  an  accurate  azimuth 
alignment. 
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After  slewing^  the  south  gyro  is  in  & favorable  position  for  accurate 
drift  determination.  The  autobiased  south  gyro  will  then  be  ready  for 
in-flight  navigation. 


3.  Off-Set  Self-Alignment 

It  is  possible  to  achieve  platform  self-alignment  with  the 
platform  coarsely  leveled  but  without  requiring  that  the  platform  axes 
be  physically  coarsely  aligned  to  north  and  east.  Instead,  the  north 
and  east,  which  are  obtained  from  an  analytic  coarse  alignment,  are 
analytically  maintained  in  the  computer.  The  fine  alignment  is  then 
achieved  by  determining  the  misalignments  between  the  computer  north 
and  east  and  the  true  north  and  east. 

Referring  to  Figure  5,  a denotes  the  off-set  of  the  platform  level 
coordinates  from  the  computer's  level  coordinates,  and  0^  is  the  azimuth 

misalignment  to  be  determined.  Under  the  off-set  condition,  the  follow- 
ing quantities  are  first  obtained: 

V„  . ad  V — Velocities  along  the  platform's  north  and  east  axes, 
NP  EP  obtained  by  integrating  the  outputs  of  north  and 
east  accelerometers 

a — The  off-set  angle,  obtained  by  a certain  coarse  alignment 
procedure,  say,  BATH. 

Next,  the  velocities  along  the  computer  north  and  east  axes  are  com- 
puted from 

V„  = V„  cos  a - v_  sin  a 
Nc  Np  EP 

(37) 

\’  - cos  a + v sin  a 

EC  EP  NP 

Finally,  a fine  alignment  technique  can  be  chosen  for  performing  the 
fine  alignment  between  the  computer  axes  and  the  earth  coordinates. 

During  the  fine  alignment,  platform  can  either  be  torqued  at  earth 
rate,  or  be  torqued  about  the  azimuth  axis  at  the  azimuth  component  of 
the  earth  rate.  The  merit  of  not  torquing  the  level  axes  is  to  avoid 
the  torquer  scale  factor  uncertainties,  which  have  been  discussed.  The 
information  of  the  torqued  coordinates  is  already  in  the  computer  since 
this  information  is  needed  to  torque  the  platform  in  the  first  place. 
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Maps! 

hhb*£«  i 
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E.  N - EARTH  EAST  AND  NORTH 
EC.  Nc  - COMPUTER  EAST  AND  NORTH 

Ep.  Np  - PLATFORM  EAST  AND  NORTH 

Figure  5*  Off-set  self-alignment. 


4.  Largs  Angle  Self-Alignment 

In  the  case  of  "large  angle  self-alignment,"  the  platform  Is 
coarsely  leveled  so  the  small  angle  approximations  for  9^  and  9^  are 

valid.  The  azimuth  axis  13  torqued  at  n^,  the  azimuth  component  of 

earth  rate.  But  the  azimuth  misalignment  9^  is  not  small  enough  to 

allow  the  use  of  small  angle  approximation. 

Under  this  condition,  a good  approximation  for  platform  drift  can 
be  obtained  from  Equations  (1)  and  (2)  as 


DS  * DN  - nAeE  + °H  “*  \ 


°E  + aAeN  “ % sin  9A 


Equations  (38)  and  (39)  can  be  solved  for  cos  0^  and  sin  0^s  respectively. 

Thus  tan  0^  can  also  be  obtained.  During  the  self-alignment  period, 

the  variation  in  azimuth  misalignment  is  small,  so  it  is  assumed  that 
©A  “ eAo*  The  result*n8  gyrocoiapassing  equation  is  therefore  given  by 

0-0  - tan”1  nAeN0  * (°E  " DE0) 

®a  ao  *:an  fl  - (T)  - D'  'l  • (40) 

A A0  *aeo  (Dn  dno^ 

The  ambiguity  of  double  values  c an  be  resolved  by  solving  0^  from 

Equation  (39)  also,  and  compare  its  sign  to  that  determined  by  Equation 
(40).  The  determination  of  0NQ,  0£O,  D^Q,  and  D^Q  will  be  discussed 

in  Sections  VI  and  VII. 

One  may  wonder  why  0^  is  not  determined  directly  from  either 

Equation  (38)  or  (39)  by  taking  the  inverse  of  cosine  or  sine  function. 
The  reason  is  that  tangent  function  possesses  steeper  slopes,  enabling 
a more  accurate  determination  of  0^. 

Notice  that  the  large  azimuth  angles  for  which  this  method  is 
intended  cannot  be  arbitrarily  large.  They  must  be  within  the  limits 
that  Equations  (1)  and  (2)  are  valid. 


5.  A Typical  Self-Alignment  Procedure 

By  combining  a few  of  the  aforementioned  concepts,  a typical 
self-alignment  procedure  can  be  developed.  As  an  example  we  may  have 
a "two-position  off-set  zero-torquing  self-alignment."  Figure  6 shows 
the  coordinates  involved  in  this  method.  In  Figure  6 N and  E are 
earth's  north  and  east;  Nc,  Ec,  and  Sc  are  the  north,  east,  and  south 

known  to  the  computer;  and  N£  and  E£  are  north  and  east  of  the  platform. 

The  90-degree  slewing  for  the  second  position  is  indicated  by  dashed 
arcs.  The  slewing  is  not  required  to  be  precise. 

The  required  alignment  steps  for  this  example  can  best  be 
described  with  the  help  of  an  activity  flow  diagram  shown  in  Figure  7. 

If  the  time  for  achieving  BATH  is  shorter  than  a coarse  alignment 
slewing,  the  total  alignment  time  can  be  reduced  when  the  off-set 
3elf-alignment  concept  is  employed.  This  is  because  the  physical 
coarse  alignment  twl.es  time  to  achieve,  while  the  off-set  technique 
90-degrees  slewing  can  be  very  rapid  without  worrying  about  its 
accuracy. 
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PLATFORM  IN  ANY  OFFSET  POSITION 


ESTIMATE  OFFSET  a BY  "BATH" 


- OPEN  LEVELING  LOOPS  AT  ACC'MTR  OU  TPUTS.  (K) 

- LEVEL  AXES  NOT  TORQUED.  w 

- AZIMUTH  AXIS  TORQUED  AT  SlA. 

- MEASURE  DELTA  VELOCITY  PULSES  FROM 
ACC'MTR  OUTPUTS. 


COMPENSATE  FOR  ALL  KNOWN  BIAS  AND  DRIFTS  (b)  j 

INTRODUCED  BY  THE  PLATFORM  AND  ELECTRONICS  FIRST 

WHICH  INCLUDE:  POSITION 

GYRO:  BIAS,  SCALE  FACTOR,  MASS  UNBALANCE.  (AUTOBIASING) 

ACC'MTR:  BIAS,  SCALE  FACTOR,  ORTHOGONALITY  I 

SYMMETRY-  I 


COMPUTE  AVg  ANDAVg  FROM  THE  COMEPNSATED  (g) 
AVn  AND  AVe  BY  COORDINATE  TRANSFORMATION 
EQUATION. 


LEAST-SQUARE  DATA  REDUCTION  TO  GET  (DJj),. 


SLEW  PLATFORM  ABOUT  90  dag. 


OBTAIN  VALUE  OF  NEW  a BY  "BA  i n.‘ 


REPEAT  FINE  ALIGNMENT  STEPS  A,  B,  AND  C. 


LEAST-SQUARE  DATA  REDUCTION  TO  GET  (0N>2. 

COMPUTE  AZIMUTH  MISALIGNMENT  USING 

(DLU-(D'n)1  t 

0 * — S-= 5LI  ♦ (0M)2  tan  L. 

A neoaL 


SECOND 

POSITION 

(ALIGNMENT) 


UPDATE  ALIGNMENT. 

TORQUE  OUT  MISALIGNMENTS. 


- READY  FOR  LAUNCH. 

- CLOSE  LEVELING  LOOPS  AT  THE 
END  OF  ALIGNMENT  PROCEDURE. 


Figure  7.  Two-position  off-set  self-alignment# 
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Section  VI.  STATE  ESTIMATION  FOR  SELF-ALIGNMENT 


1.  Platform  Alignment  State  Vector 

Platform  leveling  and  gyrocompassing  amount  to  the  determina- 
tion of  misalignments  0^,  0^,  and  0^  which,  in  turn,  require  the  know- 
ledge of  D'  and  D'.  A preferred  procedure  is  to  first  determine  0„,  9_, 
D^,  and  Dg  from  the  platform's  sensor  outputs.  The  azimuth  0^  can  then 

be  determined  using  the  gyrocompassing  Equation  (30).  In  this  procedure 
the  desired  self-alignment  state  vector  x consists  of  four  elements : 


2.  Zero-Torquing  Measurement  Equation 

Consider  the  case  of  zero-torquing  fine  alignment  where  the 
state  vector  is  determined  from  the  output  data  of  accelerometers.  For 
a sufficiently  short  fine  alignment  time  the  equations  for  misalignments. 
Equations  (27)  and  (28),  can  be  approximated  by  their  Taylor  series 
expansion  up  to  the  second  power  of  t;  i.e., 

VC)  - %0  + “nO*"  - 2 °E0  V2  <4l> 

* SE0  + “eO*  + 1 Dio  ■ <42> 


Notice  that,  since  the  north  gyro  is  not  torqued,  the  drift  term 

includes  the  north  component  of  the  earth  rate.  The  drifts  and  mis- 
alignments transform  into  accelerations  via  tilts  of  north  and  east 
accelerometers,  giving 


aN  “ -g  X sin  0£  - -g0E 
aE  * g x sin  0N  - g0N 


* 


where  g is  the  gravitational  acceleration.  Substituting  Equations  (41) 
and  (42)  into  Equations  (43)  and  (44)  gives 

®N  “ “g0EO  “ gDE0t  “ f DN0  aAfc 


®E  “ g0NO  + gDN0t  “ I °E0  nAfc  * 

Integrating  Equations  (45)  and  (46)  from  0 to  t gives  the  changes  of 
velocities  during  that  period  as 

VN(t)  " -8V  " 2 V2  * f °N0  V=3  C«) 


VE(t)  * SV  + 2 V2  ' f DE0  V3  • <48> 

Let  V « [VN,  VE]fc  be  the  measurement  vector,  the  matrix  form  of 
Equations  (47)  and  (48)  is 

M fo  -gt-fnAt3  -ft2!  LN0] 


vE(t)  gt  o 


6U  , v — A 

A 2 


ft  - f nAt3  d: 


Equations  (47)  and  (48) , or  Equation  (49)  are  the  measurement  equation 
desired. 

Grouping  Equations  (19),  (20),  (27),  and  (28)  together,  the  state 
vector  at  any  time  is  related  to  its  initial  value  by 

sin  n.t  cos  O.t  - 1 


cos  0At 

-sin  fl^t 

sin  0At 

cos  fi.t 

Equations  (49)  and  (50)  provide  a way  for  determining  the  desired  state 
vector  from  the  integrated  measurement  data. 


It  may  be  asked  Why  the  state  vector  is  not  determined  using 
Equations  (45)  and  (46)  as  measurement  equations.  The  answer  is  that 
the  integration  reduces  errors  caused  by  the  quantization  effect  of  the 
accelerometer  output. 


In  practice,  discrete  measurements  are  made  at  time  t * xk  for 
k » 1 ~ N,  and  x is  the  sampling  period.  Substituting  the  discrete 
time  into  Equations  (49)  and  (50)  gives 


Vk>" 

0 

-gxk 

-fvv 

- f tV  ‘ 

vE(W 
» ■ 

g'*k 

0 

K 2,  2 
“ T k 
2 K 

- f n.-rV 

D A 

sin  0. xk  cos  n.xk  - li 


Vk>l 


1 - cos  n.xk  sin  n.xk 


Vk> 


D'(lc) 


0 0 C08  n.Tk 

A 


•sin  nATk 

A 


D'OO 


0 0 sin  nAtk 


C08  n.Tk 
A 


If  the  total  measurement  time  is  sufficiently  small. 


sin  0Axk  « 0Axk 


cos  n.xk  - 1 
A 


then  an  approximation  for  Equation  (52)  is 

feN(k)l  fl  0 xk  0 


*E<k> 


n* 

N(k) 


-vki 


D*(k) 


o o n.xk 

A 


sss 


3.  Estimation  Technique* 

When  the  environment  is  ideal , where  disturbance  and  noise 

are  not  present,  measurement  of  V (t)  and  V_(t)  at  two  different  values 

N E 

of  t is  sufficient  for  determining  the  state  vector.  In  reality,  the 
measurements  are  contaminated  by  noise  due  to  ground  vibration,  wind 
buffeting,  instrument  noise,  and  other  random  disturbances.  Under  this 
condition,  accurate  determination  of  the  state  vector  demands  the  use 
of  a large  number  of  redundant  measurements  in  conjunction  with  a 
statistical  estimation  technique. 

Two  well  known  approaches  of  statistical  estimation  techniques 
applicable  to  platform  self-alignment  are  the  least  square  regression 
approach  and  the  Kalman  filtering  approach.  A comparison  of  the  two 
approaches  is  given  as  follows: 

a)  Least  Square  Regression  - This  approach  does  not  make  use  of 
statistical  properties  of  the  noise,  if  available.  The  associated  data 
reduction  process  can  be  made  either  sequential  or  batch.  If  the 
desired  number  of  state  estimations  Is  much  less  than  the  number  of 
measurements,  the  least  square  regression  algorithm  requires  less  com- 
puter time  and  smaller  computer  memory  as  compared  to  Kalman  filtering 
14]. 

b)  Kalman  Filtering  [5]  - The  Kalman  filtering  algorithm  has  a 
built-in  provision  for  taking  advantage  of  known  second  order  statistics. 
The  associated  data  reduction  process  is  sequential,  which  generates  an 
estimate  of  the  state  from  each  measurement. 


We  choose  th ' least  square  regression  approach  for  the  platform 
self-alignment.  The  choice  is  based  on  two  facts:  first,  knowledge 

of  the  statistics  of  the  noise  is  not  good  enough  to  enjoy  the  merit  of 
Kalman  filtering;  and  secondly,  the  required  number  of  state  estimations 
is  far  less  than  the  number  of  measurements  so  least  square  regression 
approach  is  superior  in  the  required  computer  time  and  memory. 

In  Section  VII  a new  least  square  regression  algorithm,  tailored  to 
our  platform  self-alignment  application,  is  presented  in  detail. 


26 


Section  VII.  A NEW  LEAST  SQUARE  ALGORITHM 
1.  Measurement  Equations 

Least  square  regression  method  requires  a special  form  for 
measurement  equations.  Rearranging  Equation  (51)  and  adding  measurement 
noise  to  it,  we  get 

2 3 

\ . t !_*■  » a /i_\  /rev 


Where  n^k)  and  ng(k)  are  additive  noise,  k * 1 ~ N,  and 

*i  - 

A a .&  D* 

A2  2 °E0 

*H-fDH0V3 


+ Vk) 

(55) 

+ nE<k> 

(56) 

* 1 ~ N,  and 

A3  ' *8N0T 

*4  " f "NO'2  • <58> 

\ -6  °E0  V3 

The  problem  becomes  the  determination  of  A^,  Ag,  A^,  and  A^  from  a 
large  set  of  redundant  measurements  VN(k)  and  V£(k).  and  A^.  are  not 
needed  because  they  differ  from  A^  and  A^,  respectively,  only  by  a 
known  constant  multiplier. 


2.  Tha  Usual  Least  Square  Algorithm  [67,8] 

Measurement  Equations  (55)  and  (56)  are  In  the  form  of  three- 
term  third-order  polynomials  in  k.  A conventional  set  of  least  square 
regression  formulas  are  available  in  textbooks  in  statistical  mathematics 
for  determining  the  coefficients.  Applying  conventional  formulas  to 
our  problem,  the  coefficients  are  determined  from 
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A2 
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Y2 

Y3 

• 
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A3 

-1 

Z1 

A4 

- C 

Z2 

\ 

.Z3 

N 


'j-2 


kJ  VN(k) 


j - 1 ~ 3 


k-1 

N 


Zj  “ 2 ^ VE<k)  J “ 1 ~ 3 


k-1 


N 


:i  - 1 


k-1 


(59) 


(60) 


(61) 


(62) 


(63) 


(64) 


These  formulas  have  been  used  for  platform  self-alignment  as  well 
as  their  applications  and  are  applicable  to  any  three-term  third-order 
polynomial.  The  algorithm  can  be  written  either  for  batch  processing 
or  for  sequential  processing. 


3.  A New  Least  Square  Algorithm 

Examining  Equations  (55)  through  (58),  it  is  seen  that  they 
can  be  expressed  in  the  following  equivalent  form: 

VN(k)  - Aj>  + A2k2  + uA4k3  + n^(k)  (65) 
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Vg(k)  - A3k  + A^k  - uA2kJ  + ng(k)  (66) 


where 


(67) 


Notice  that  some  of  the  coefficients  are  correlated  deterministically. 
Can  this  property  be  employed  to  improve  the  estimation  accuracy?  The 
answer  is  in  the  affirmative. 


Let  us  study  a more  general  case 

VR  - Ak  + Bk2  + uDk3  + Ug 

V£  - Ck  + Dk2  + vBk3  + n£ 


(68) 


where  VR  and  are  measurements  made  at  sampling  instants  k ■ 1 ~ N, 

and  A,  B,  C,  and  D are  parameters  to  be  estimated.  A cost  function  is 
chosen  as  follows: 


(69) 


The  cost  function  is  to  be  minimized  by  an  optimum  selection  of  Ay  B, 
C,  and  D.  Setting 


Si  . Si  Si  Si 

Ia  “ Sb  “ Sc  3d 


o 


(70) 


a set  of  four  algebraic  equations  are  obtained  which  can  be  combined  into 
a single  matrix  equation: 


where 


" C2 

C3 

0 uC, 

4 

C3  C4 

' uC£ 

1 

g 

■p* 

o 

G « 

0 

-uC4 

C2  C3 

uC. 

4 

0 

C3  C4  + uC6 

N 

Ci“ 

i - 2,  3,  4,  and  6 

k-1 

N 

N 

H1  - I kVN 

(k). 

w * 

w2 

S k2[VN(k)  - ukV  (k)] 

k-1 

k-1 

N 

H 

w.  - \ kvE(k), 

W s 
w4 

£ k2[VE(k)  + ukVjj(k)  1 

k-1 

k-1 

Values  of  A^,  i ■ 1 ~ 4,  are  given  by 


> 

1 

— 1 

A2 

-1 

W2 

» G 

A3 

W3 

i 

1 

.V 

(77) 


(78) 


(79) 


(80) 


Notice  that  the  measurements  V„(k)  and  V_(k)  contribute  to  the  estima- 

tion  via  the  computation  of  W^,  i - 1 ~ 4.  Also,  the  property  of  the 

platform  kinematics  together  with  the  choice  of  cost  function  result 
in  the  absence  of  C<.  in  the  computation. 


The  estimate  for  the  initial  state  vector  is  obtained  from  A., 

x 

i - 1 ~ 4,  using  Equations  (57)  and  (58). 


(82) 


Consider  a simple  case  having  a single  measurement  equation 
VK  - Ak  + Bk2  + Bk3  + . 

Let  A ■ 4,  B - 2,  and  k * 1 ~ 10.  The  values  of  VR  are  generated  by 

adding  noise  n^  to  the  value  of  polynomial  at  each  k.  The  values  of  n^ 

are  taken  from  a table  of  normally  distributed  random  numbers,  having 
zero  mean  and  a variance  of  one.  These  values  are  listed  in  Table  1. 

TABLE  1.  MEASUREMENT  DATA  FOR  EXAMPLE  1 


K 

\ 

1 

1 

6.724 

2 

31.628 

3 

-1.377 

82.623 

4 

2.334 

178.334 

5 

318.864 

6 

528.414 

7 

-0.494 

811.506 

8 

1.048 

1185.048 

9 

0.347 

1656.347 

mm 

0.637 

2240.637 

Both  new  and  usual  least  square  algorithms  are  used  to  estimate  A 
and  B from  the  data  VR,  k ■ 1 ~ 10.  The  results  are  shown  in  Table  2, 

listing  values  of  estimates  and  their  percentage  error  as  compared  to 
true  values.  It  is  apparent  that  the  uew  algorithm  produces  much  more 


TABLE  2.  RESULTS  FOR  EXAMPLE  1 


A « 4 


B - 2 


Usual  Algorithm 
N » 10 

New  Algorithm 
N - 10 

New  Algorithm 
N - 5,  (1st  5) 

New  Algorithm 
N « 5,  (2nd  5) 


3.67969 

2.09863 

8% 

5% 

3.95581 

2.00111 

1% 

0.05% 

3.87604 

2.00359 

3.1% 

0.17% 

4.00342 

2.00059 

0.086% 

0.025% 

accurate  estimates.  The  new  algorithm  gives  better  results  even  if 
fewer  measurements  are  used.  Appendix  A contains  the  computer  program 
for  this  example. 

Example  2 

Consider  the  following  case  of  two  measurement  equations: 

V « Ak  + Bk2  + Dk3  + tl.  \ 

2 3 } • (83) 

« Ck  + Dk  + BkJ  + n£  I 

Let  A ■ 4,  B * 2,  C ■ 3,  D ■ 1,  and  k * 1 ~ 10.  The  values  of  and 
Vjj  are  generated  In  a manner  similar  to  that  for  Example  l.  Table  3 
lists  the  measurement  data. 

TABLE  3.  MEASUREMENT  DATA  FOR  EXAMPLE  2 


”k 

°K 

VK 

V* 

K 

-1.276 

-1.218 

5.724 

4.782 

-0.318 

-0.799 

23.682 

25.201 

-1.377 

-1.257 

55.623 

70.743 

2.334 

-0.337 

114.334 

155.663 

-1.136 

0.642 

193.864 

290.642 

0.414 

-0.011 

312.414 

485.989 

-0.494 

0.364 

468.506 

756.364 

1.048 

0.037 

673.048 

1112.034 

0.347 

2.816 

927.347 

1568.816 

0.637 

0.563 

1240.6’? 

2130.j63 

Again,  both  new  and  usual  least  square  algorithms  are  employed 
to  estimate  A,  B,  C,  and  D.  The  results  are  listed  in  Table  4,  with 
the  percentage  error  of  each  estimate  estimated.  Again,  the  new 
algorithm  gives  much  better  results.  Appendix  B presents  the  computer 
program  for  this  example. 

TABLE  4.  RESULTS  FOR  EXAMPLE  2 


Usual 

Algorithm 

New 

Algorithm 


3.61133 

9.7% 

3.95093 

1.2% 


2.08984 

4.4% 

2.00272 

0.14% 


2.07031 

31.0% 

2.88883 

3.7% 


1.00099 

0.1% 


4.  Sensitivity  Consideration 

It  is  expected  that  the  new  algorithm  is  less  sensitive,  as 
compared  to  the  usual  algorithm,  with  respect  to  erratic  measurement 
data  and  to  computation  errors.  The  rationale  is  that  the  coupled 
parameters  have  a tendency  to  hold  each  other  at  their  nominal  values. 
Example  3 will  show  this  effect. 

Example  3 

This  example  uses  the  same  measurement  model,  same  measurement 
data,  and  same  computer  programs  as  those  used  for  Example  2.  To 
observe  the  effect  of  erratic  measurement  data  on  estimates,  the 
measurement  V(10)  is  changed  by  1%  and  a set  of  new  estimate  for  A, 

B,  C,  and  D is  made  using  new  and  usual  algorithms.  To  observe  the 
effect  of  computation  error  on  estimates,  the  value  of  V(10)  is  restored 
to  its  original  value  and  the  value  of  Cg  is  changed  by  0.01%.  Another 

set  of  estimates  are  made  using  both  algorithms.  All  estimates  are 
listed  in  Table  5.  Values  of  sensitivity  in  the  table  are  calculated 
using  the  formula 

„ New  estimate  - Nominal  estimate 
* Nominal  estimate  * ' ' 

The  result  confirms  the  expectation  that  new  algorithm  has  lower 
sensitivity. 


In  this  example  a real  world  platform  system  is  considered.  The 
platform  is  of  the  class  proposed  for  PERSHING  II  application.  Twelve 
hundred  and  fifty  pairs  of  measurement  data  were  recorded  at  outputs 
of  north  and  east  accelerometers.  The  total  sampling  time  is  240  seconds 


TABLE  5.  SENSITIVITY  COMPARISON  FOR  EXAMPLE  3 


Condition 


Algorithm 

Used 


Nominal 

Usual 

New 

3.61133 

3.95093 

2.08984 

2.00272 

2.07031 

2.88883 

1.22070 

1.00099 

17.  Change 
in  V(10) 


0.017.  Change 
inC6 


Usual 


49.8% 


34.97. 


-14.47.  -26.77. 


-0.057. 


-1.37. 


15.47. 


-93.67. 


52.57. 


-0.067. 


2.877. 


-0.057. 


A - Estimate  of  A 
S - Sensitivity 

The  quantities  to  be  estimated  are  misalignments  and  platform  drifts. 
New  and  usual  least  square  algorithms  are  used  for  data  reduction.  The 
former  consists  of  Equations  (77)  through  (81)  while  the  latter  consists 
of  Equations  (59)  through  (64)  and  (87).  The  results  of  data  reduction 
are  shown  in  Table  6.  To  explore  the  sensitivity  of  both  algorithms 
with  respect  to  computation  errors,  a poor  matrix  inversion  subroutine 
is  used.  When  computations  are  done  with  double  precision,  both 
algorithms  produce  reasonable  results.  When  computations  are  done  wii-h 
ordinary  precision,  the  result  from  new  algorithm  is  not  reasonable, 
knowing  the  quality  of  the  platform  used.  But  the  result  from  the 
usual  algorithm  is  ridiculous  regardless  of  the  platform  considered. 

The  results  show  the  superiority  of  the  new  algorithm. 


TABLE  6.  RESULT  FOR  EXAMPLE  4 


Ordinary  Precision 

Double  Precision 

New  Algorithm 

Usual  Algorithm 

New  Algorithm 

Usual  Algorithm 

1162.8  arcsec 

4.736  X 10^  arcsec 

-10.55  arcsec 

-6.294  arcsec 

695.6  arcsec 

2.599  X 10^  arcsec 

691.6  arcsec 

693.0  arcsec 

0.275  deg/hr 

5.362  X 1036  deg/hr 

13.306  deg/hr 

13.19  deg/hr 

-0.0017  deg/hr 

2.824  X 1036  deg/hr 

0.0775  deg/hr 

0.038  deg/hr 

The  only  disappointment  in  this  example  is  that  the  exact  values 
of  misalignments  and  drifts  were  not  available  at  the  time  of  experiment, 
therefore  a precise  comparison  of  two  results  could  not  be  made.  To 
partially  overcome  this  difficulty,  theoretical  error  analyses  are 
developed  in  Section  VIII.  These  analyses  will  help  to  evaluate  the 
quality  of  algorithms.  Appendix  C presents  two  computer  programs  used 
in  this  example. 
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Section  VIII.  THEORETICAL  ERROR  ANALYSES 


The  approach  of  theoretical  error  analyses  used  here  is  to  develop 
analytical  relationships  relating  the  standard  deviation  of  estimation 
error  to  the  standard  deviation  of  the  noise.  The  analyses  are  done  for 
new  and  usual  least  square  algorithms. 


1 . Analysis  for  New  Algorith  m 

Recall  Equation  (80)  and  define 


H » G 


H is  a 4x4  matrix  whose  ij-element  will  be  denoted  by  h^ . Expanding 

the  first  row  of  Equation  (86)  and  using  the  relationships  in  Equation 
(79) , Aj  can  be  expressed  as 

.4 

= 1 wj 

j=l 

■ hn  I kvk>  +>’i2  !l  kVk>  - » | k3vk>! 

k U k I 

+ h13  £ kVE(k)  + + u Y k3VN(k)|  . (87) 

k Ik  k / 

The  error  in  A^  is  caused  by  errors  in  V^(k)  and  V^(k)  which  are  denoted 
by  e^(k)  and  e^k),  respectively.  Let  e^  represent  the  error  of  A^, 
then  from  Equation  (87)  we  can  get 
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h “ hll  1 keN(k)  + h12(  1 k\W  - U S k3eE(k)i 
k k k J 

+ k13  S l“B(k)  +hu(S  k2'E(k)  +u  £ k3eH<k)| 
k ( k k I 

■ S |(huk  + ht2k2  * h14ak5)eH(k) 

+ (hi3k  + hHk2  -hn“k3K<k)} 


+ ^h13k  + h14k“  - h12uk-'jeE(k) j (88) 

which  expresses  the  error  in  in  terms  of  source  errors. 

Assume  the  following  statistical  properties  for  source  errors: 

a)  Zero  mean,  i.e., 

< eN(k)  > - < eE(k)  > - 0 (89) 

where  "<  >'  denotes  ensemble  average 

b)  Uncorrelated  between  axes  and  from  time  to  time,  i.e., 

< eN(i)eN(j)  > - < eE(i)eE(j)  > - 0 i*j  (90) 

< eN(i)eE(j)  > - 0 all  i,j  (91) 

c)  Equal  and  stationary  variance,  i.e., 

< e2(k)  > » < e2(k)  > » <r2,  u constant.  (92) 

N Hi  c 

Taking  the  ensemble  average  over  the  square  of  Equation  (88),  applying 
previously  mentioned  error  properties,  and  rearranging  terms,  we  can 
obtain  the  error  variance  of  as  follows: 


°L  -°e  2 |(hllk  + h12k2  + uh14k3) 

k > 2| 

+ (h13k  + hl4k2  - uh12k3) 


We  shall  digress  for  a moment  to  derive  a number  of  relationships 
which  will  help  to  simplify  the  final  expression.  Substituting  details 


eatsfe  uoauima  t&aas>o& < 


of  V„(k)  and  V_  (k)  as  given  in  Equations  (65)  and  (66)  into  Equation 
N E 

(67)  and  regrouping  terms. 


*1  - 41  1 k(hlllt  + h12k2  + Uh141'3) 


+ A2  S >2(hnk  + hi2k2  + Uhi4l‘3)  " uk3(hnk  + hi4k2  • uhi2k3j 

k 


I k(h13 


k + h^k  - uh12l< 


+ A4  1 k2(h13k  + h14k2  ' 0h12k3)  + “k3(huk  + h12k2  + uk14k3) 


comparing  both  sides  of  Equation  (94)  shows  that  the  coefficient  of 

should  be  1,  and  those  of  A2,  A3,  and  A^  should  be  zero.  Therefore  we 
get  the  relationships: 

\ 

I k(»Uk  + h12k2  + uh14k3)  - 1 
k " 

£ jk2(hllk  + h12K2  + ^l'/)  * Uk3(h13k  + hl4k2  “ uh12k3))  = 0 


1 k(hi3k + v2  - uhi2k3)  - 0 r5) 

k 

£ I kZ(h13k  + h34k2  ' al,12k3)  k uk3(hllk  + h12k2  4 l,h14k3)}  = 0 • 

k * 

Return,  to  our  analysis  of  error  variance  and  apply  the  relation- 
ships of  Equation  (95)  to  Equation  (93),  The  result  is  a very  neat 
expression, 

r.{  - -r  - . <9« 


where  is  the  normalized  error  variance  for  A^.  The  normalization  is 
done  with  respect  to  the  variance  of  source  error » 


In  the  similar  manner, 
and  are  obtained  as 


the  normalized  error  variances  for  A^,  A^, 


2 

^2 


(97) 


(98) 


(99) 


Notice  that  values  of  h^,  i * 1 ~ 4,  depend  solely  on  N,  the 

number  of  measurements,  and  u,  the  correlation  parameter.  Therefore, 

2 

normalized  error  variances  N « 1 ~ 4,  are  independent  of  measure- 
ment  data,  but  depend  on  the  kinematics  of  the  platform  misalignment. 


By  taking  the  square  roots  of  Equations  (96)  through  (99), 
equations  for  normalized  standard  deviation  of  estimates  are  obtained 
as 


1i  mJST  » 1 “ 1 ~ * » (100) 

another  very  neat  expression. 

Equation  (100)  is  very  useful  in  several  ways.  For  a given  set 
of  error  standard  deviations  of  the  source  and  the  number  of  measure- 
ments, it  can  be  used  to  estimate  the  error  standard  deviations  of 
estimates.  Hcwover,  for  a given  set  of  source  error  standard  deviations 
and  a set  of  prescribed  standard  deviations  for  estimates,  it  can  be 
used  to  determine  the  minimum  number  of  measurements  needed.  Finally, 
knowing  the  error  standard  deviation  of  the  estimate  and  the  number  of 
measurements,  the  equation  can  be  used  to  determine  the  standard  devia- 
tion of  source  error,  a tool  for  identification. 


2.  Analysis  for  Usual  Algorithm 

Recall  Equations  (59)  and  (60),  and  define 


i 


Q is  * 4 X 4 matrix  whose  ij-element  will  be  denoted  by  q^.  Adopting 

an  approach  of  derivation  similar  to  that  for  the  new  algorithm,  the 
expression  of  normalized  error  variances  for  estimates  of  A.,  i « 1 ~ 4, 
can  be  obtained  as  1 


"*i  “ qii  » 1 " 1 ~ 4 * 


(102) 


Similarly,  the  normalized  standard  deviation  expression  is  given  by 


\ qii  » 1 = 1 ~ 4 


(103) 


In  this  case,  depends  only  on  N,  the  number  of  measurements,  but  not 
on  the  correlation  parameter. 

3.  Comparison 

For  a given  measurement  condition  the  accuracy  of  estimates 
produced  by  new  and  usual  algorithms  can  be  compared  by  comparing  their 
standard  deviations.  Direct  comparison  of  Equation  (100)  and  Equation 
(103)  in  their  literal  forms  is  difficult.  A numerical  example  will  be 
used  to  demonstrate  the  superiority  of  new  algorithm. 

Example  5 

Consider  the  same  platform  alignment  problem  of  Example  4.  Figures 
8 and  9 show  plots  of  normalized  standard  deviations  as  functions  of  N, 
the  number  of  measurement.  Axes  of  the  plots  are  in  log-scale. 

For  N **  1250,  the  new  algorithm  gives 


3 G * \ 

NO  E0 


0.24  X 10 


0.26  v 10 


(104) 


while  the  usual  algorithm  gives 


n “ V,  « 0.15  * 10 
"NO  E0 

'D1  * qD'  * 0.94  X 10 

NO  E0 


(105) 


Comparing  Equation  (104)  to  Equation  (105),  the  superiority  of  new 
algorithm  is  evident. 


Appendix  D contains  the  computer  program  used  for  providing  plotting 
data  for  Figures  8 and  9, 
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NORMALIZED  ERROR  STANDARD  DEVIATION 


Section  IX.  RECOMMENDED  FURTHER  STUDY 


There  are  several  problems  which  deserve  further  study.  Solution 
to  these  problems  will  allow  a truly  optimum  implementation  of  the  IMU 
self-alignment  systems. 


Statistical  theory  shows  that  more  accurate  estimates  are  obtained 
with  larger  N,  the  number  of  measurements.  However,  computer  error 
analysis  shows  that  larger  N results  in  more  computer  error  because 
more  computation  is  involved.  It  is  desirable  to  choose  an  N such  that 
the  total  error  is  at  its  minimum*  A method  for  making  such  a choice 
is  yet  to  be  developed. 


Even  though  the  new  least  square  algorithm  is  less  sensitive  to 
computation  errors  as  compared  to  the  usual  algorithm,  it  is  still 
desirable  to  keep  computation  errors  as  small  as  possible,  expecially 
those  occurring  during  matrix  inversion.  Notice  that  the  G matrix  of 
Equation  (77)  is  symmetric  and  has  four  zero  elements.  This  special 
form  may  allow  the  development  of  a matrix  inversion  subroutine  which 
is  more  efficient  in  computation  accuracy  and  computation  time. 

The  new  algorithm  reported  here  is  given  in  the  form  of  batch 
process.  This  algorithm  can  be  modified  to  become  a sequential  process 
or  a hybrid  process  which  is  a semi -batch-semi-sequential  process. 

It  is  desirable  to  verify  the  analytically  predicted  superiority 
of  the  new  least  square  algorithm  for  platform  alignment  by  a precision 
hardware  IMU  which  can  be  calibrated  for  experimental  comparison. 

The  analytic  results  obtained  from  this  study  provide  insights 
for  coarse  alignment  which  is  required  prior  to  the  fine  alignment. 

The  coarse  alignment  can  also  be  performed  automatically  and  rapidly 
with  the  aid  of  the  computer  already  available  for  fine  alignment. 

It  will  be  interesting  to  explore  the  possibility  of  a combined  coarse 
and  fine  alignment  using  the  same  equipment  and  giving  a best  overall 
alignment  result. 
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Section  X.  CONCLUSIONS 


An  original  contribution  o£  this  study  is  the  development  and 
analysis  of  a new  least  square  regression  algorithm  specially  for  the 
self-alignment  of  IMU  systems.  It  was  shown  experimentally,  as  well 
as  analytically,  that  this  new  algorithm  is  superior  to  the  usual 
algorithm  in  accuracy  and  in  sensitivity. 

Although  the  new  algorithm  was  originally  intended  for  the  self- 
alignment of  a gimbaled  platform,  it  can  be  used  for  the  alignment  of 
a strapdown  platform  as  well,  with  some  minor  modifications.  It  can 
also  be  used  for  an  IMU  consisting  of  electro-optical  sensors,  because 
the  underlying  kinematic  principle  is  similar. 

Other  results  of  this  study  include  thorough  derivation  for  drift 
equations,  misalignment  equations,  and  the  gyrocompassing  equation. 
Several  self-alignment  concepts  were  reviewed  and  discussed  using  the 
analytic  foundation  developed.  Five  examples  were  developed  to  help 
in  confirming  the  theoretical  prediction. 

Several  areas  deserving  further  investigation  were  recommended. 
The  solution  to  these  areas  will  allow  a truly  optimum  implementation 
of  IMU  self-alignment  systems. 


45 


3 


REFERENCES 


Hung,  J.  C.,  State-of-the-Art  Self-Alignment  Techniques  for  Fixed 
Base  Inertial  Platforms.  Report  for  Task  Order  72-625,  US  Amy 
Missile  Comnand,  Redstone  Arsenal,  Alabama,  December  1973, 

White,  H,  V.,  Serco  Analysis  of  an  Inertial  Platform,  Report  No. 
RG-TR-63-3,  US  Army  Missile  Command,  Redstone  Arsenal,  Alabama, 
July  1963. 

Pitman,  G«,  Editor,  Inertial  Guidance.  John  Wiley  and  Sons,  New 
York,  1962. 

Hung,  J.  C.,  A Study  of  Several  Estimation  Techniques  for  Gyro- 
compassing.  A Report  for  Task  Order  74-237,  US  Army  Missile 
Command,  Redstone  Arsenal,  Alabama,  May  1974. 

Kalman,  R.  E«,  "A  New  Approach  to  Linear  Filtering  and  Prediction 
Problems,"  Trans.  ASME,  J.  Basic  Eng..  Vol.  82D,  1960,  pp.  34-45. 

Hildebrand,  F.  B.,  Introduction  to  Numerical  Analysis.  McGraw-Hill 
New  York,  1956. 

Ralston,  A.,  A First  Course  in  Numerical  Analysis.  McGraw-Hill, 

New  York,  1965. 

Burr,  I.  W.,  Applied  Statistical  Methods.  Academic  Press, 

New  York,  1974. 

H aiming,  R.  W.,  Numerical  Methods  for  Scientists  and  Engineers, 
McGraw-Hill,  New  York,  1962. 

i.  Cannon,  R.  H.,  "Alignment  of  Inertial  Guidance  Systems  by  Gyro- 
compassing  — Linear  Theory,"  Journal  of  the  Aerospace  Science. 
Vol.  28,  No.  11,  1961,  pp.  885-895  and  912. 

Kayton.  M.,  and  Fried,  W.  R.,  Avionic  Navigation  Systems, 

John  Wiley  and  Sons,  New  York,  1969. 

Danik,  B.,  Gyrocompass  Alignment  of  Inertial  Platform.  Document 
No.  KD-73-9,  The  Singer  Co,,  Kearfott  Division,  6 February  1973. 


^igp^^lPp!il||iP|S§P55R|SBESBliE 


Appendix  A.  COMPUTER  PROGRAMS  FOR  EXAMPLE  1 

(IN  BASIC) 


10  READ  C2>C3,C4,C5,C6 
20  DATA  385^3025,25333.220825 . . 

40  LET  fll*C2 

41  LET  Q2*C3+C4 

42  LET  Q3-C4+2*C5+C6 
50  LET  0*01*03-02*02 

60  LET  Ml =Q3/D 

61  LET  M2=-Q2/D 

62  LET  M3=M2 

63  LET  M4*Ql/D 
100  DIM  VC101 

110  FOR  2*1  TO  10 
120  READ  VCU 
130  PRINT  "V("I")»"VC13 
140  NEXT  I 

150  DATA  6*724,31 •6288,82*623, 178*334,318.864 

160  DATA  528.414,811*506,1185*05,1656.35,2240.64 

180  LET  XI  =0 

190  LET  X2=0 

200  FOR  J=l  TO  10 

210  LET  Xl=J*VCJ3+Xl 

220  LET  X2*J*J*<l+j)*VC*n+X2 

230  NEXT  J 

260  LET  A=M1*X1+M3*X2 
270  LET  3=M3*X1+M4*X2 
280  PRINT  "A=”A,"3="3 

290  PRINT  "TRUE  VALUES  ARE:  A«4  3*2  C*B*2” 

300  END 
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10  READ  -2*C3*C4*C5*C6 

20  DATA  385-3225,25333,220605., 1 T97RaaR+a* 

105  REM  COMPUTATION  OF  K1  TO  K9 

110  LE»'  D«G2*CC4*C6-C5*C5>-C3*CC3*C6-C4*C5)+C4:*’CC3*C5-C4*C4> 

140  LET  K1=<C4*C6-C5*C5)/D 

150  LET  K2*CC4*C5-C3*C6)/D 

160  LET  K3*<C3*C5-C4*C4)/D 

170  LET  K4*K2 

180  LET  K5*(C2*C6-C4*C4)/D 

190  LET  K6»(C3*C4-C2*C5>/D 

200  LET  K7*K3 

210  LET.  K8*K6 

220  LET  K9“<C2*C4-C3*C3)/D 

300  REM  READ  IN  OBSERVATIONS 

305  DIM  VC103 

310  FOR  1*1  TO  10 

320  READ  VCI1 

330  PRINT  ,,V<MI*,)*MVCI  3 

340  NEXT  I 

358  DATA  6.724,31  .682*82. 623.-  ‘78*334,318. 864 

360  DATA  528.414,811.506,1185.05*1656.35*2240.64 

400  REM  WITH  PRECOMPUTED  CONSTANTS*  DATA  PROCESSING  BEGINS  HERE 

405  LET  XI  *0 

406  LET  X2=0 

407  LET  X3-0 

410  FOR  J*l  TO  10 
420  LET  X1*U*VCJ]+X1 
430  LET  X2«J*J*VC J1+X2 
440  LET  X3*J*U*U*VC«J1+X3 
458  NEXT  J 

500  LET  A«K1*X1+X2*X2+K3*X3 
510  LET  B*K4*X1+K5*X2+K6*X3 
520  LET  C*K7*Xl ♦K8*X2+K9*X3 
530  PRINT  ”A*,,A*,,B*,,B*”C*,,C 

550  PRINT  ’'TRUE  VALUES  ARE:  A*4  B*2  C*B=2,, 

600  END 


Appendix  B.  COMPUTER  PROGRAMS  FOR  EXAMPLES  2 AND  3 

(IN  BASIC) 


10  PRINT  "LEAST  Svi.  ALGO.  USING  PARAMETER  CORRELATION' 

12  LET  02=385 

13  LET  03=3025 
U LET  04=25333 

15  LET  C5=220625. 

16  LET  C6=l .97&40E+06 
18  PRINT 

20  DIM  VC103 
30  FOR  1=1  TO  10 
40  READ  VC  I 3 
62  NEXT  I 

70  DATA  5.724,23.682,55.623, 1 14.334, 193.864 

72  DATA  312. 414,466. 506,673. 048,927. 347,1240. 64 

80  DIM  UC10  3 

93  FOR  J=1  TO  10 

100  READ  UCJ3 

120  NEXT  J 

130  DATA  4.782,25.201,70.743* 155. 6c:, 290. 642 
132  DATA  485.989,756.364,1112.04,1568.82,2130.56 
200  DIM  UC4, 1 3 

210  LET  VC1, I ]=VC2, 1 3=VC3, 1 3=WC4, 1 3=0 

240  FOR  K= 1 TO  10 

250  LET  WC1,1 3«WCl,l 3+K*VCX3 

260  LET  V C 2 , 1 3 = V C 2 , 1 3 +K*K*  < V C X 3 +K*U C K 3 > 

273  LET  U C 3 , 1 3 = V C 3 , 1 3 +X*J C K 3 

2 80  LET  WC4,  1 3=v.’C4,  1 3 +K*K* CUCK3 +K*V CK3 ) 

293  N^.-.T 

430  DIM  G C4, 4 3 

401  LET  GCl, 1 3=GE3,33=C2 

432  LET  GC1,23=GC2, 1 3=GC3,43=GC4,33=C3 

433  LEV  GC 1 ,33=GC3, 1 3=0 

404  LET  GCl, 43=GC4, 1 3 =G C2 , 3 3=G C3, 23=04 

405  LET  GC2,2 3=G C4, 4 3 =C4+C6 

406  LET  GC2,43=GC4,23=2*C5 
460  DIM  H C 4 , 4 3 

470  HAT  H=INV(G) 

5C0  DIM  XC4, 1 1 

505  PRINT  "X  VECTOR" 

506  PRINT 

510  MAT  H=H*V.’ 

520  MAT  PRINT  X 

5 to  PRINT  "TRUE  VALUES:  XI =4  X2=2  X3=3  X4  = l" 


-vr^— r.v  r S'  .r,  ■u^u*-^*^***^^ 


[?£* , : - "•' 

^g?;,jrT--ig 

p£ffi?3b^-'Ca 


PRINT  "THE  USUAL  LEAST  S3.  ALGO.” 

DIM  VC103 
FOR  1=1  TO  10 
HEAD  VCI3 
NEXT  I 

5,724>23. 682*55. 623*  i 14.334*  193.864 
DIMAUcl0J4l/4'468*506'673‘348>927'3/‘7',243#64 

FOH  J=1  TO  10 
HEAD  UCU3 
NEXT  J 

DATA  4.782*25.201*70.743* 155.663*293.642 
DATA  485.989*756.364* 1 1 1 2 .04* 1 568 .82*2 1 33 .56 
HEAD  C2*C3*C4*C5*C6 

DATA  385*3025*25333*220825.* 1 .97840E+06 

LET  D=C2*<C4*C6-C5*C5)-C3*(C3*C6-C4*C5>+C4*(C3*C5-C4*C4> 

LET  K1=<C4*C6-C5*C5)/D 

LET  K2=K4=<C4*C5-C3*C6)/D 

LET  X3=K7=<C3*C5-C4*C4)/D 

LET  X5=<C2*C6-C4*C4)/D 

LET  K6=X8-(C3*C4-C2*C5) /D 

LET  K9=<C2*C4-C3*C3)/D 

LET  Xl=X2=::3=0 

FOH  1=1  TO  10 

LET  X1=X1  + I*VCU 

LET  X2=X2+I*I*VCI3 

LET  X3=X3+I*I*I*VCI 3 

NEXT  I 

LET  YI=Y2=Y3=0 
FOH  u=l  TO  10 
LET  Yl»Yl+J*UCJ3 
LET  Y2=Y2+u*«j*UCJ3 
LET  Y3=Y3+J*J*«J*UCJ3 
NEXT  J 

LET  A=Kl *X1 +K2*X2+K3*X3 
LET  3=K4*X1+X5*X2+K6*X3 
LET  C=H7*Hl  ♦i\8*X2+K9*X3 
LET  E=X1*Y1+K2*Y2+K3*Y3 
LET  F=K4*Y1 +X5*Y2+K6*Y3 
LET  G=K7*Y1+X8*Y2+K9*Y3 
PH I NT  ”A="A*"3="3*"C="C 
PH  I NT  "E=”E*  "F="F* "G="G 
PH  I NT 

PHIKT  "THUE  VALUES  AHEs" 

PH  I NT  "A=-  ”*”S=2”*"C  = r* 

PHINT  ••E=3”*••F  = iM*,•G=2,, 
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Appendix  C.  COMPUTER  PROGRAMS  FOR  EXAMPLE  3 

(IN  FORTRAN) 


PROGRAM  MA IN ( INPUT .OUTPi IT «TAPER=INPUT ,TAPF6*0UTFUT> 

C**«**LEAST  SQUARE  ALGORITHM  I iS I NO  PARAMETFR  CORRFt  AT I ON 

INTEGER  OELPM»  OELPE 

OIMFNSION  C(S)»  G(4,4)»  0ELPN(1?50)»  OFLPE(1250).  VN(1?S0)» 

1 VFUP50).  W<4),  X (A)  « VNOFF  { 1250)  » VEOFF(1250>,  H(4,4> 

DOUBLE  PRECISION  G,  H.  W,  X 

C***«*SETTIMG  UP  C(2>  TO  C(6)  < C(R>  IS  NOT  USED,  ) 

DO  ? 1=2,6 
C<I>  = 0.0 
00  4 K=1 , IPSO 
FK=K 

4 C(I)=C(I).FK**I 
2 CONTINUE 

WPITE  <6,6)  (I,  C ( T » * I=?«6> 

6 FORMAT  ( 1H1  ///  S(10X.*C*n*=*E20,12/)) 

C***«*FSTAPLISHING  G-MATRIX 

C 
C 
C 
C 
C 
C 


FARTH  PATE  = 7.2921 1F-05  RAO I AN/SF CONO 
LATTTUOP  = 34.6425  DEGREES 
SAMPLING  PERIOD.  TdlJ  = 0.1R2  SECOND 
U = -(7-COMPONENT  OF  EARTH  RATE) *TAU/3 

= - (-7.2921  IE-05  « SIN34.6425)  * 0.19?  / 3 
= 0.?CS29477P?E-6S 
U=0.?6S?9473«2E-0S 
G 1 1 » 1 ) =C (?) 

G ( 1 ,2) =C (3) 

G(1.3) =0.0 
G ( 1 ,4) =U*C (4) 

G (?, 1 ) =G ( 1 »?) 

G ( 2 , ? > =C ( 4 ) -U*C ( 6 ) 

G (?«  3)  = (-U) *0(4) 

G(?.4)=0.0 

G(3,l)=0«0 

G(3«?)=G(?»3) 

G( 3.3) =C (?) 

G 3.4)=C(3) 

G(4,l)=G(l»A) 

G(4,?) =0.0 
G (4,3) =G  < 3.4) 

G (4,4 ) =C (4) *U*C(6) 

WRITF  ' P>  ( (G ( T I)  , J=1.4>»  1 = 1,4) 

« FORMAT  ( \H0,7X, *0-MATPI X»/4X  »4F?4. 14/4X  »4E?4, 14/ 
l 4X,4P?4. j 4/4X.4F24. 14 ) 


r a»Bo#riFTTING  G = G- INVERSE  (IN-PLACE  STORAGE) 

CA Ll  MTXINV(G,4) 

WRITF(6,9)  ((G(J,J)»  J=l,4).  1=1*4) 

9 FORMAT  (1H0/7X,*G-TNVERSF*/4X,4F24.14/4X,4E24.14/ 
\ 4X,4F?4,\4/4X,4F?4.14) 
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C-1} 


C»»**«»»»*«FOR  DFRtJGGING  USE  ONLY 

C TO  CHECK  THE  TNVEP5F  OF  THE  INVERSE  OF  G 

DO  *0  1=1*4 
00  «1  *1=1  *4 

81  H(I*J)=G(I*J) 

80  CONTIMJF 

CALL  ' TXINV(H*4) 

*‘RITE(6.82)  ((H(I,.i>.  J=l*4),  1 = 1*4) 

82  FORMAT  (1H0/7X**G-IniVF»SE',S  INVERSF*/4(4X.4E24.14/) ) 

C**********DEBUGGIMG  INFORMATION  ENOS 

C*****PEAO  IN  MEASURED  OAT* 

REA0(5» 10)  (DELPNU).  DflpE(I).  1=1*1250) 

10  FOPMAT (10(214)) 

WRITE (6*1?)  (DELPN(I) » 1=1,1250) 

12  F0PMAT(1H1*10X**0ELPN1I)  « 1 = 1 TO  i?50*t//63 (8X,2()I6/) ) 
WPITE<6,13)  (DELPE ( I ) « T = 1*12S0) 
n FORMAT (1H1*10X**0ELPE cl ) • 1=1  TO  1250**//63<8X*20I6/) ) 

C*****COMPUTING  VNOFF(t)  AND  wFOFF ( I ) (IN  NUMREP  OF  PULSES) 
VNOFF ( 1 ) =OELPN U ) *55 
VEOFF(l)=OFl.PE(l)*40 

no  40  1=2*12^0 
AUXM  = OFLPN(I)  ♦ 55 
AUXE  = DFLPF(I)  ♦ 40 
VNOFF ( I ) =VNOFF U-l ) *AUXN 

40  veoff(I)=.vfoff<i-ii*auxf 

WRITE (5*42) 

42  FORMAT  (1H1,10X*«I**8X*«VNOFF (I)»*12X,*VEOFF(I)4) 

WRITE' (6  *44)  (25*  I.  VNOFF  (?5*  I)  « VE0FF(25*D*  1 = 1*50) 

44  FORMAT (4X. IB, 2E24. 14) 


1 c*****OBT4INING  VNU)  A NO  V£  U ) R 1 COOPOIANTF  TRANSFORMATION 

* C UN  NUMBER  OF  PULSfS) 

C HATH  « 150  n^GRFES  = 2.5)7993878  PAOIANS 

f PATH=2.fcl 7993878 

| SR=STN (HATH) 

g C9=CO$(HATHJ 

= DO  14  1=1,1250 

I VN ( I ) = VNOFr(I)*rB  - VF0FF(I)*SR 

: 14  VE(I)  = VFOFc-(I)»CH  ♦ VMOFFU)*SR 

I WRITE (6*30) 

I 30  FOPMAT  UHl,10X,*I*,IOX.*VNU>»*15X,*VFU)*) 

| WRITE {4*15)  (25*1.  VN(25*T).  VF(?5*I),  1=1*50) 

fc  IS  FOPMAT  (4X.TR.2F 24 *14) 


C***««C0MP'ITING  W-VECTOP 
W(1)=0.0 
W (2) =0.0 
W(3)=0.ft 
W(4)=0.0 
DO  16  K = 1 . 1 2S  0 


52 


FK=K 
W(  1 ) 
W(?) 
W<3> 
16  W<4) 


«M1)  = W<1)  ♦ FK*VN(K> 

w (?)  = w(?>  ♦ FK*r«f*(VN(K»-u«FK*ye<Kn 

W<3>  = W(l>  ♦ FK*VE(K> 

W(4)  s Wf 4)  ♦ FK*FK*( VE (K) ♦U*FK*VM(K) ) 

WPITE (6» 17)  (I*  HMD.  I = l*4> 

FORMAT  ( 1 HI  ///  4 (]4Xf *W#I1*=*E?0,1?/)  ) 


C*«  »«*COMPl  IT  I NO  X-VECTOP 
00  18  1=1 »4 
> ( I) - C . 0 
no  ?0  J=1 *4 

PO  X < I > = X(I)  ♦ 0 ( I • J) *w ( J) 

1H  CONTINUE 

WPIT E(6*??)  <I«  X ( I ) • 1=1*4) 

;>?  FORMAT  ( 1 HO //  4 ( 14X*»X*T! *=*E?0, l?/> > 

C«^***COMPlJTING  PLATFORM  PARAmfjERS 
C N-AXIS  SCALE  FACTOd  SFK=1 00441, 

C F-AXIS  CCALE  FACTOR  $FE=101712. 

SFN=100441, 

SFE=10171?. 

TAU=0.19? 

Arl./TAU 
B=?./TAU**P 
ZETANO=A*X (3) /SEN 
7FTAF.0=  (-A)  *X  ( 1 ) /SF£ 

ORFTMO-H*X (4) /SEN 

ORFTFO=(~R)*X(?)/SFR 

WRITE  (6«?4)  7ETANQ*  ZtTaFO,  ORFTNO,  ORFTEO 


?4  FORMAT  UHO//)4X,*7FTANn=*E?0,lP* 
\ i4x,-»7F.TAEn=*F?0, 1?* 

1 14X,fiORFTNn=:*E?0.1?* 

1 14X,  ^r)RFTEn=«E?0.1?* 

C<f**««ENO  OF  THE  ESTIMATION  PPOC-PAm 
END 


PARIAN*// 

RADIAN*// 

PAOIAN/SEC*// 

RADIAN/SFC*) 


igl  spi ' 
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SUBROUTINE  MTXINV(G.M) 


ON  OUTPUT  6 IS  THE  INVERSE  OE  6 
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i~- 

10 

I?;- 

20 

§L 

E 

s 

! 30 

1 

ke 

40 

50 

1 

60 

f 

70 

• 

80 

l 

90 

i 

i 

100 

no 

tf 

k\r- 

120 

| 

130 

1 

140 

C«««««0N  INPUT  G IS  G« 

DIMENSION  G(4,4) 

DOUBLE  PRECISION  G 
DO  140  K=1,M 
IF  <G(* *K) ) 10*  160* 

GZZ=1 .Q/G(K  »K) 

DO  90  I— 1 .M 

IF  (I-K)  20.  90*  60 
CONST=G(I,K)*GZZ 
DO  SO  J*I*M 
IF  (J-K)  SO*  50.  40 
G( I , J) =G( I. J)  ♦CONST*G  < J*K) 

GO  TO  50 

G(I.J)=G(I.J)-CONST*G<K,J) 

CONTINUE 
GO  TO  90 

CONST *G (K  * I ) *GZZ 
00  BO  J=T.M 
IF  (J-K)  170,  80.  70 
G ( I ♦ J) =G ( I . J) -CONST*G (K  * J) 

CONTINUE 
CONTINUE 
DO  UO  J=K.M 
IF(K-J)  100.  110.  180 
G(K« J) =G(K, J)*GZ7 
CONTINUE 
00  130  1=1. K 
IF  (K-I)  190,  130.  120 
G ( I *K) = (-G( I • K ) ) *G77 
CONTINUE 
G(K,K)=G7Z 
CONTINUE 
DO  ISO  I=?,M 
J1 =1-1 

00  ISO  J=!,J1 
G ( I . J) =G ( J. I ) 

PE TURN 
WRITE (S.P10) 

GO  TO  POO 
WPITE(S,??0) 

GO  TO  ?00 
WRITE (6.P30) 

GO  TO  POO 
wPITE(*,?40) 

CONTINUE 

FORMAT (4X,*7PP0  01  AGON Al  ELEMENT.  BAD 
FORMAT (4X»*EPR0R  IN  INDEXING  IN  DOING 
FORMAT (4X,*f PROP 
FORMAT <4X,*ERR0R 
END 


ISO 


160 


170 


180 


190 

200 

P10 

2P0 

?30 

240 


IN 

IN 


INDFXING 

INDEXING 


IN 

IN 


DOING 

DOING 


DATA*) 

DO  BO*) 
DO  110*) 
DO  ,30*) 
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PROGRAM  MA  IN  ( INPUT .OUTPl IT . TAPES* INPUT « T APE6»0UTPUT > 

C THE  USUAL  LEAST  SQUARE  ALGORITHM 
C 3-TEPM  3RD  ORDER  POLYNOMIAL 

INTEGER  DELON*  OELPE 

DIMENSION  OELPNt 12*50) ♦ nFLPE(l?50)»  VN(1250)«  VE(1250>» 

1 VNOEE (1250) ♦ VEOEE (1250) • 0(3.3) 

DOUBLE  PRECISION  O 

C*****SETTING  UP  C2  TO  Cfi 
C2=0«0 
C3=0.0 
C4=o.n 
cs=o.o 
C6s n.o 

DO  2 K=l.l?50 
FK=K 

C2=C?*EK**? 

C3=C3'*FK**3 
C4=C4*FK*»4 
C5=C5*FK**S 
? C f,=Cfi*FK**G 

WRITF(G*4)  CP*  C3.  C4*  CS.  C6 
4 FORMAT  (1HI///10X.*C2=*f?0.12/10X.*C3=*E20.!2/ 

I 10X,*C4=*F20.12/IOX.*C5=*F20.12/10X.»C6=*E20.12) 

C**««*COMPUTING  D AND  01  TO  OR 
0(1*11 =C2 
0(1 *2) =C3 
0(1*3) =C4 
0(2*1) =C3 
0(2.2)=C4 
0(2*3) =CS 
0(3*1) =04 
0(3*2)=CS 
0(3*3) =Cfi 
CALL  SYMINV(0*3) 

WRITE  (6*6)  ( (Q(I«.J)  « J=)*3)*  1 = 1*3) 
f>  FORMAT(1HO/7X*«0-INVE»SF*/3(4X.3E24.14/) ) 

C«*«»opEAO  IN  MFASIIREO  DATA 

READ (S* 10)  (DELPN(I) « OFLPE(I>«  1=1*1250) 

10  FORMAT (10(214)) 

WRITF (ft, 1?)  (OELPN(I)*  1=1*1250) 

1?  FORMAT ( 1H1 * 10X**0ELPN( I ) « 1=1  TO  1250*. //63 (BX. 2016/) > 

WRITE (4* 13)  (DELPE(I) » 1=1*1250) 

13  FORMAT(1H1,10X.*OELPE(I) . 1=1  TO  1250*»//63(flX. 2016/1) 


C«ofl«»C()MPUTING  VMOFF(I)  ANO  VFOFF(I)  (IN  NUMBER  OF  PULSES) 
VNOFF  ( l ) =0El.pN  ( 1 ) *55.0 
VEOFF ( 1 ) =0FLP£ ( 1 ) *40*0 
DO  40  I=2»1?S0 
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VNOFF  < I ) sVNOFF ( 1-1) *OELPN ( t ) *55.0 
40  VEOFF ( I ) =VEOFF ( 1-1 ) ♦OELPF ( l > *40.0 

4?  F0PM4T<S ( IHl , 1 OX  * *1 • « 8X  « *VNOFF  ( I > * * \2* ** V' E Ui FF ( I > *> 
WR1TF (4*44)  (25*1 « VNOFF (25*1) . VE0FF(25«I>*  1-1»S0) 

44  F0RMAT(4X«I8«2E24.14) 


•OBTAINING  VN ( I ) ANO  VE C I * BY  COOROIANTE  TRANSFORMATION 

(in  number  of  pulsfsi 

RATH  = 150  OEGRFES  * 2.G1799387B  RADIANS 
BATH=?.617RB387R 


Sf?=SIN (BATH) 

CB=COS(BATH) 

no  14  1=1*1250 

VN(  T)  = VNOFF ( I } *CB  - VEOFF(I>«SR 

14  VE(I)  = VEOFF (I)*CB  ♦ VNOFF ( 1 ) *SB 

50  FORMAT  (lHl,10X.*I*.lOX,«VN(I)».15X.*VF(n*) 
WRITF(fi,15)  (25*1*  VN(25*I)»  VF(25*1).  1=1*50) 

15  FORMAT  (4X«IB*2E24.14) 

COMPUTING  VI*  Y2*  Y3*  2)*  22*  23 

Y1=0.0 

Y2=0.0 

Y3=0.0 

00  17  1=1*1250 
FI=I 

Y1=Y1*FI*VN(I> 

Y2=Y2*FI*FI*VN(I> 

17  Y3=Y3*FI*FJ*rl*VN( I) 

71=0.0 

22=0.0 

Z3=0.0 

00  1«  J=l»1250 
FJ=J 

Zl=71 *FJ*VF ( J) 

7?s7?*FJ*F J*VE ( J) 

IB  73=Z3*FJ*FJ*FJ*VF(.J) 


••♦•COMPUTING  XI  TO  X4 

X 1 =0 ( l « M *Y1 *Q ( 1 *2) *Y2*o ()*3) *Y3 
X2=0(2« l )*Yl *Q (2*2) *Y2*0(?*3) *V3 
X3=0( 1 *1) *7 1*0 <1*21 *Z2*0(1 »3>*73 
X4sO(2. 1 )*Z1*0(2*2' *22*0 t 2*3) *73 
VR ITF (G*PO)  «1*X2*X3»X4 

20  FORMAT  (IHI///1AX.*X1**F20.12/HX.*X?=*E20.1?/ 
I 14X**X3=*P20.12/1AX«*X4=*E20.12) 


•COMPUTING  PLATFORM  PARAMETERS 

M-AXI5  SCALE  FACTOR  SFN*100441, 
E-AXIS  SCALE  FACTOR  SFE*10lT12. 
SFN= 100441. 

SFE=101712. 

TAU=0.1R? 


! 


A=1./TAU 
P=2./TAU**2 
ZETAN0=A*X3/SFN 
ZETAE0=<-A)*X1/SFE 
0RFTN0=8*X4/SFN 
0PFTE0=(-R)*X2/SFE 

WRITE  (*>«?4)  7ETAN0.  ZETAFO*  DRFTNO*  ORFTEO 
24  FORMAT  < 1H0//14X«*ZETANOv*E20. 1 2*  RADIAN*// 

1 14X.*ZETAEOr»F20.12*  RAOTAN*// 

1 14X«*ORFTNn=*E20.12*  RAOI AN/SEC*// 

1 14X«*ORFTEO=*E20.12*  RADIAN/SEC*) 

j C*****END  OF  THF  FORTRAN  PROGRAM 

END 
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SUBROUTINE  SYMINV(G.M) 


ON  OUTPUT  6 IS  THE  INVERSE  OF  G 


C*****ON  INPUT  G IS  G, 

DIMENSION  G(3,3) 

OOURLE  PRECISION  G 
DO  140  K=1,M 
IF  (G(K,K) ) 10*  160.  10 
10  3ZZ=1 »0/G (K.K) 

DO  90  1=1. M 
IF  (I-K)  20.  90.  60 
20  CONST*G(T,K)*GZZ 
DO  SO  Jal.M 
IF  (J-K)  30.  50.  40 
30  G(I«J)=6(1«J) *CONST*G ( J.K) 

GO  TO  SO 

40  G(I.J)=G(I.J)-CONST*G<K,J) 

50  CONTINUE 
GO  TO  90 

60  CONST =G (K , I ) *G2Z 
DO  80  J=I.M 
IF  (J-K)  170.  80.  70 
70  G(I.J)=G(I.J)-CONST*G<K,J) 

80  CONTINUE 
90  CONTINUE 

DO  110  J=K.M 
IF(K-J)  100.  HO.  180 
100  G(K,J)=G(K,J)*GZZ 
110  CONTINUE 

00  130  1=1, K 
IF  (K-I ) 190,  130,  120 
120  G(I,K)=(-G<I,K))*GZZ 
130  CONTINUE 
G(K,K)=G7? 

140  CONTINUE 

DO  ISO  1=2, M 
J1=I-1 

DO  ISO  J=1,J1 
150  G(I* J)=G( J.I) 

RETURN 

160  WRITE (6,210) 

GO  TO  200 
170  WRITE (6.220) 

GO  TO  200 
180  WRITE (6.230) 

GO  TO  200 
190  WRITE (6,240) 

200  CONTINUE 

210  FORMAT (4X,*7ER0  DIAGONAL  ELEMENT,  RAD 
220  FORMAT (4X.*EDR0P  IN  INOFXING  IN  DOING 
230  FORMAT (4X,*EOROR  IN  INDEXING  IN  DOING 
240  FORMAT (4X,*FRR0P  IN  INDEXING  IN  DOING 
END 


DATA*) 

DO  80*) 
DO  110*) 
DO  130*) 
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Appendix  0.  COMPUTER  PROGRAMS  FOR  EXAMPLE  5 

(IN  FORTRAN) 


PROGRAM  MA I N { I NPUT .OUTPl IT . TAPES* INPUT .T APE6*0UTPUT ) 
C*#o*#TARLF  OF  NORMALIZED  STANDARD  DEVIATIONS  FOR  THE  USUAL 
C AND  THE  NFW  LEAST  SQUARE  ESTIMATION  ALGORITHMS 

o I MENS I ON  A (A.4) • 011(30).  022(30)*  033(30)* 

1 0(4.4).  Mil (30) . H22(30) * H33 (30) « H44(30) 

•1=0  . ?f>G?R473a2E»0S 
C?=0.0 
C3=0.0 
C4=ft.O 
C5=0.0 
C6=G.O 

no  ?0  N=l«30 
M= (N-l ) *100*1 
I=N*1 00 
no  30  j-m.t 
c?=c?*j**?.o 

C3=C3*J**3.0 
C4=C4*J**4.0 
CS=CS*J**S.O 
30  C6=C4*j«*iS.O 

A ( 1 . 1 ) =0? 

A(1.?)=C3 
A ( 1 .3)  CA 
A (?« 1 ) =03 
A(?.?)=C4 

A(2.3)=C* 

A(3.D=C4 
A ( 3 ♦ P ) =CS 
A(1.3)=C* 

CALL  MTXTNV(A.l) 

Oil (N)=A(l«l) 

0?2(N) =A (?.?) 

033 (N) =A (3.3) 

G(1.1)=C? 

G( 1 *?) =C3 
G(l.3)=0.0 
G(1«4)=U»C4 
G (2. I ) =G  < 1 ♦?) 

G(Z.?)=C4-II«C6 
G(?*3)=(-U)*C4 
r,(?,4)=0.0 
G (3. 1 ) =0.0 
G(3«?)=G(?.3) 

G(3.3)=C? 

G(3«4)=C3 
G(4.1)=GC1 *4) 
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G<4.?)=0.0 

G(4.3)=G(3,4) 

G(4*4)=C4-*U*C6 
CALI.  «TXTMV(G.4) 

HU  fM)sA(l«n 

H33<N)=GC»*3) 

?0  H44  (N!)  =G  (4.4) 

•«VITP(^.40) 

40  FOPM4T  < 1*1/1  0X.»U/100*.19***CU  <N)  *. ’4*  ,#Hl  1 (N)  *.?5X  .*H33  (N)  *) 
WPITF(ft,90)  (W,  nll(\*).  HI  1(N>*  H31(W),  N=  i * 30  > 

^0  POP'*AT  ( lfiK.IG.  10X«C?0.  IP.IO^.PPO.1?*  1 OX.c?0.12) 

•rfPITF  ( ^ * A,  0 > 

AO  FORMAT  ( 1H1  /10X.«-N»/100»«19x.»C??(Ni)*.?4x,*H??(N)**?SX**H44(N)*) 
wPTTF(A.70)  (N.  02?(N)«  H??(N).  H44(N).  N=1.30) 

7-  FORMAT  ( ]0X,lc,i0X*F?0*l2.10<*F?n.l2*l0X.F2(l.l2) 

•*,CIITF  (A. 30) 

40  FORMAT  ( 1*1  /]  OX* *M/1  00*»1C)X**C3'MN)'B') 

'•'PITF(A.qO)  (.'J*  M=l*30) 

q?  FOOMAKlOX.T^.lOX.PPO.l?) 

FMP 
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IwPUT  G IS  6.  ON  OUTPUT  G IS  THF  T‘ IVFRSE  OF  G 
DIFFUSION  G(4«4) 

no  140  k=1«m 
IF  (G(K,K))  10*  160*  10 
10  GZ7=1.0/G(K*tO 
00  90  1=1 ,M 
IF  (I-K)  20,  90.  SO 
20  CONST=G<T.K)*GZ7 

no  so  j=i.m 

IF  (J-K)  30.  50.  40 
10  0(1*  J)=G(l.J)  ♦CONST*G ( J,*) 

00  TO  SO 

40  0 ( I , J ) =0 ( I . J) -COMST*G IK,  J) 

SO  CONTINUE 
00  TO  QO 

SO  C0MST=O(K.T) »OZ7 
00  40  J=  I « ’■* 

IF  (J-K)  170.  00.  70 
70  0 ( I , J)  =0  ( I , J)  — CONST *0  (K  • J) 

«0  FONT  INI  IE 
90  CONTINUE 

DO  110  J=K , M 

ifo<-j)  ion,  no.  iso 
100  0(K,J)=0(K,J)*GZ7 
110  CONTINUE 

no  no  1=1,  k 

IF  (K-I)  190.  no.  120 
no  G(I,K)  = (-G(I.K)  )*G7Z 
110  CONTINUE 
G(K,K)=G7Z 
140  CONTINUE 

OO  ISO  I=?,M 
J1=I-1 

no  iso  j=i,ji 

ISO  G(I,J)=G(J.I) 

RETURN 

ISO  WRITE IS, 210) 

GO  TO  200 
170  WHITF(S,??0) 

GU  TO  200 
nO  WRITF(S.?10) 
no  TO  200 
190  '.'WITF(S,240) 

200  CONTINUE 

2)0  FORMAT (4X,»7FR0  DIAGONAL  FI.E^FNT • SAD  nATA*) 

2?0  FORMAT (4X,*FRR0P  IN  INOFXlNG  IN  DOING  DO  90*) 

210  FORMAT (4X,*FRR0p  TN  INOFXlNG  IN  DOING  00  l 10*) 

R40  FORMAT (4X,*FPR0R  TN  INOFXlNG  IN  DOING  DO  130*) 

End 
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